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Abstract 

We present a computer program which solves the Schrodinger equation of the 
stationary states for an average nuclear potential of Woods-Saxon type. In this 
work, we take specifically into account triaxial (i.e. ellipsoidal ) nuclear surfaces. 
The deformation is specified by the usual Bohr parameters. The calculations are 
carried out in two stages. In the first, one calculates the representative matrix of 
the Hamiltonian in the cartesian oscillator basis. In the second stage one diagonalizes 
this matrix with the help of subroutines of the EISPACK library. If it is wished, one 
can calculate all eigenvalues, or only the part of the eigenvalues that are contained 
in a fixed interval defined in advance. In this latter case the eigenvectors are given 
conjointly. The program is very rapid, and the run-time is mainly used for the 
diagonalization. Thus, it is possible to use a significant number of the basis states 
in order to insure a best convergence of the results. 
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Licensing provisions: none 

Computer: PC. AMD Athlon 1000MHz 

Hard disk: 40 Go 

Ram: 256 Mo 

Swap file: 4 Go 

Operating system: WINDOWS XP 

Software used: Microsoft Visual FORTRAN 5. OA (with full optimizations in 
the settings project options) 

Programming language: fortran 77/90 (double precision) 

Number of bits in a word: 32 

Number of lines: 3150 lines with comments 

Keywords: Nuclear physics, Energy levels, Wave functions, Schrodinger equa- 
tion, Woods-Saxon potential 

Nature of the problem: 

The Single particle energies and the single particle wave functions are calcu- 
lated from one-body Hamiltonian including a central field of Woods-Saxon 
type, a spin-orbit interaction, and the Coulomb potential for the protons. 
We consider only ellipsoidal (triaxial) shapes. The deformation of the nuclear 
shape is fixed by the usual Bohr parameters 7) . 

Method of solution: 

The representative matrix of the Hamiltonian is built by means of the Carte- 
sian basis of the anisotropic harmonic oscillator, and then diagonalized by a 
set of subroutines of the EISPACK library. 

Two quadrature methods of Gauss are employed to calculate respectively the 
integrals of the matrix elements of the Hamiltonian, and the integral defining 
the Coulomb potential 

Restrictions: 

There are two restrictions for the code: 

The number of the major shells of the basis does not have to exceed Nmax=26. 
For the largest values of Nmax (~23-26), the diagonalization takes the major 
part of the running time, but the global run-time remains reasonable. 
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Typical running time: 

(With full optimization in the project settings of the Microsoft Visual Fortran 
5. OA on Windows XP ) 

With NMAX=23, for the neutrons case, and for both parities, if we need all 
eigenenergies and all eigenfunctions of the bound states, the running time is 
about 80 sec on the AMD ATHLON computer at 1GHz. In this case, the cal- 
culation of the matrix elements takes only about 20 sec. 
If all unbound states are required, the runtime becomes larger. 

Long write-up 



1 Purpose of the Fortran program 



1.1 The Schrodinger equation 



The program solves the Schrodinger equation for one body-deformed potential: 

= Etlh) (1) 

Here H represents the Hamiltonian of the system (neutrons or protons), and, 
Ei , and, 0j , represent respectively, its eigenenergies, and its eigenfunctions. 



1.2 The Hamiltonian 



The Hamiltonian of the nucleon is defined by [1,2]: 

H = T + V + V so + e<f) c (2) 

The quantities T, V, V s0 , indicate respectively, the kinetic, potential , and 
spin-orbit energy. For the proton, the Coulomb potential is represented by 
C , and e is its charge 
Explicitely: 

T = -f.V (3) 

2m v ; 

h = Planck constant 
m= nucleon masse 

Owing to the fact that the nuclear forces have a short-range character, the 
average nuclear potential must " follow on average" the nuclear density dis- 
tribution: 
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For the case of the spherical symmetry (see also (19)) , we have: 



V{r) [<x - ; i£W = 1 , °,r-R, W 



l + exp(^)J l + eX p(^) 
For the deformed case, the above definition is generalized as [9] : 



Vn 

V(r) = — (5) 

K J 1 + exp(R v Lyjay) v ; 

with 

Vo, Ry, fly = mean field parameters 
Ly = quasi-radius (see eq. (9)) 

V so {r) = ~ (VS(r) Ap)a (6) 

with 

p = (fi/i) V = Neutron or proton momentum 
a = (cr x ,o-y,cr z ) = Pauli spin-matrices 

For the same reasons that for V, the mean field in the expression of the V so 
operator is given by a similar definition (see also the subsection 3.1) : 

S(f) = " - - — - (7) 

1 + exp(i? so L so /a so ) 

with 

k = spin-orbit coupling strength (there, the quantity is absorbed by k, this 
latter is integrated to 5'(T > ) and expressed in MeV.fm 2 ) 
R so , a so = mean field parameters of S(f) 

L so is the quasi-radius of the spin-orbit mean field (see eq. (12)) 



2 The Coulomb potential: 



For the protons, the Coulomb's potential is approximated by the one of the 
liquid drop model [1,2] 



$ C (Z,P,$) = 



Pcharge 



Z 2 



dzx 



Zi 



2tt 



z -Z) d -§ + 2p 2 s - 2 Ps Pcos(ip - <£>) - 2P^ sm(ip - $) 



dip 



yj(z - Zf + p 2 s + P 2 - 2p s P cos(v? - $) 



(8) 
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where: 

(Z, P, $) = cylindrical coordinates of the point where the Coulomb potential 
is calculated, (here, Z must not be confused with the protons number) 
Pcharge = — l)e/(4/3)7i\R^= charge density of the liquid drop 
(Z — 1)= "number of protons in the liquid drop" 
R ch = radius of the charge density 

The integration domain is defined by the volume limited by the surface n c h = 
(see eq. (15)). 

The "nuclear surface of the protons" is given by p s = Psurface m ^ ne ec L- (16). 
In fact, the code computes directly the quantity e$ c (which is the Coulomb 
energy of the proton in the Coulomb field) instead the Coulomb potential $ c .( 
see the function ephi. in the code). 



3 Supplementary details on the deformation of the mean field, and 
the different nuclear surfaces: 



3. 1 The quasi-radius and the nuclear surfaces 



In the central average potential, the information on the distortion of the nu- 
clear surface is given by the dimensionless quasi-radius Ly(f). which is defined 
as [1]: 

MO = - n Jf5> 0) 



Ri 



vn y (r) 



The quantity Hy(r) is defined so that to recover the well-known spherical case 
(see section 4). 

riy(r) = \jlt V (f) ~ TTVmin - y/-^Vmm (10) 

Here, 7Ty m i n is the absolute minimum of ny(r). This latter describes an hyper- 
surface which is not the real nuclear surface, because generally 7Ty(^) 7^ in 
the expression of the quasi radius. The actual nuclear surface may be obtained 
by putting ix v (~r) = 0. 

In this work, we have restricted ourselves only to simple ellipsoidal (triaxial) 
shapes for the effective nuclear surface. 

Mr1 = |r + |r + ^-i = o (ii) 

Ay, By, and, Cy are thus the semi-axes of the ellipsoid. 
In fact, we have to consider three distinct interactions, i.e. the central, the 
spin-orbit, and the Coulomb interaction. Therefore, we must define three re- 
spective surfaces (equations (11), (14), (15)). Thus, the equation (11) describes 
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the nuclear surface relatively to the central interaction V(r). 
In completely analogous way we have to define similar quantities for the spin- 
orbit interaction 



L so {r] 



vnjr 



(12) 



with, 

n«o( r ) = \J ftsoij) ~ ft so mm ~ V — ft so mm (13) 

Where, 7r somin is the absolute minimum of 7r so (r), and the effective "spin-orbit 
surface" is written as: 

«so(r) = ^r + j^ + -^-l = (14) 

^-so n so °so 

For the Coulomb potential, the effective nuclear surface is defined in the same 
way: 

2 2 2 

ft ch (?) = 4r + |2- + 7^-i = o (15) 

^-ch n ch ^ch 

Following the expression of the equation (8), the Coulomb potential must be 
expressed in cylindrical coordinates Therefore, the equation of the effective 
"Coulomb nuclear surface" (15) can be rewritten as: 



fl 2 = 1 ~ {Z/C ch f 

r surface f, / A \ 2 i / • I T> \ 2 \ V ' 

{{cos ip/A ch ) + (sin ip/Bch) j 



where: 



x = Psurface COS (p , y = p surface Sin if , Z = Z 

The "three densities" (neutrons, protons, spin-orbit) differ very little from each 
other. Therefore, the three surfaces are homothetic and slightly different to 
each other. Nevertheless, in order to simplify the problem, the protons distri- 
bution is assumed to be uniform in the calculation of the Coulomb potential. 



3.2 The deformation parameters 



Since we have three similar surfaces, and, so as to avoid repeating three times 
the same thing, we will omit to specify the indices of the surfaces. For example, 
the three volume conservation conditions are simply replaced by only one 
equation : 

^nABC = ^nR 3 (17) 

Actually, because of this condition, only two parameters are necessary 

As usual, one will prefer the Bohr parameters ((3, 7) instead of those of the 
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ellipsoid. 



B 



C 



R 

X 
R 

X 
R 

X 



/ c x 1/2 



V4W cos ^-r ] 



cos(7 - -tt) 



/ 5 x 1/2 

/ 5 \ 1/2 
! + /?(-) cos( 7 ) 



(18a) 
(18b) 
(18c) 



i? is the radius, and, x insures the volume conservation condition (17), 



4 The Spherical case. The mean field parameters 



4-1 The case of spherical symmetry 



When A=B=C , or when the Bohr parameter (3 = 0, the nuclear surface 
becomes spherical, and Ry, R so , or, R c h represents simply the nucleus radius. 
In this case, we obtain the familiar Fermi-function for the two mean potentials 
((5),(7)): 



V(r) 



1 + cxp [(r - R v ) /a v ] 



S(r) 



K 



1 + exp [(r - R so ) I a s 



i.e. potentials of Woods-Saxon type. 

The Coulomb's potential (8) reduces to the well-known form: 

$ c (r) = [(Z - l)e/2R ch ] [3 - (r/R ch ) 2 ] if r < R ch 



$ c (r) = {Z - l)e/r if r > R, 



ch 



(19) 



(20) 



The spin-orbit interaction (6) can be expressed in the spherical case as: 

V°° = -(^1 T. A = A tt (21) 



dr r r dr 

Finally, the V so operator takes the familiar form, 

r or 



(22) 



The relations (19-22) of the spherical case are not used explicitly in the code. 
However, the well known spherical degeneracy of the energy levels were used 
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in order to check the program. Moreover, the relation (20) serves as a first 
checking for the Coulomb potential. 



4-2 Mean field parameters 

Two options have been included in the code in order to choose between a 
particular set of parameters, or the Myers parameters [3]. Thus, it is possible 
to define its own parameters in a separate file, or to employ those of Myers. In 
this latter case, the calculations are made in a suitable subroutine. In fact, only 
a part of the parameters set, namely, V , Ry, R so , Rch , is deduced from the 
droplet model of Myers, the remaining, i.e. k, ay,a so , are extracted from the 
Ref. [1]. The explicit expressions for these parameters are given in appendix. 



5 Principle of resolution 

5.1 The method 

The principle of this method consists to look for the eigenfunctions of the 
Schrodinger equation by their expansion on a truncated basis of the harmonic 
oscillator. In other words, the method used in solving such problem amounts 
essentially to writing the representative matrix of the Hamiltonian in this 
basis. 

In practice, this method is characterized by two distinct stages. First, one 
builds the representative matrix of the Hamiltonian by means of the cited 
basis. Next, one diagonalizes this matrix in order to obtain the eigenvalues 
and the eigenvectors. 

In our work, the cartesian coordinates are the most suitable. 



5.2 The harmonic oscillator basis 

The basis functions of the harmonic oscillator are defined as: 

\n x n y n z T) = i ny <p nx (x)4 ny (y)4 riz (z).a^ (23) 

There, i ny is a phase factor which insures, in accordance with the imposed 
symmetries (see section 6), that the matrix elements are real. 
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Explicitly: 



^ ( x ) = V A* ex P [- (/W 2 / 2 ] • ^ ( A^) Ac = 

with analogous expressions for the y, and z axes. 
The intrinsic spin states are: 



(24) 



0+1/2 = 



0-1/2 = 



(25) 



The /i na: , (or /i„ , or /i n J quantities symbolize the normalized Hermite poly- 
nomials. 



hnM = H nx (x)/yJ(2"*.n x \irV*) (26) 
H Hx (x) are the usual Hermite polynomials 

The quantum numbers integers, and give the order of the 

Hermite polynomials, E = ±| represents the projection of the intrinsic spin 
on the z-axis. 

At last, m and u x , u y , cj z represents the mass of the oscillator, i.e. the mass 
of the nucleon, and, its frequencies. 

For convenient, the quantities ftw x , frw y , fkv z ,(or respectively /3 x ,/3 y ,f3 z ) are 
called the deformation parameters of the basis. 

If the three frequencies are equal, the oscillator is then isotropic, and it can 
be characterized by only one frequency (cu x = u y = u z = u ) 
Furthermore, we assume that the nuclear surface of the oscillator is an equipo- 
tential. This involves the following condition: 

uj x uj y uj z = cJq or {fv-jJx) ■ (twy) ■ (twz) — (^o) 3 (27) 



5. 3 The representative matrix of the Hamiltonian 



With help of this basis, the matrix elements of the Hamiltonian H can be 
written as: 



H \n x n y n z T) 



n' x n' y n' z X> 



T + V + V so + e.$ c \n x n y n z T) (28) 
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c 



5.3.1 Matrix elements of the mean field V and the Coulomb energy e.$ 

Since V does not depend on the spin, we adopt the following convenient no- 
tation for V: 



V \n x n y n z T) = (n' x n' y n' z \V\ n x n y n^j x <5 S 's (29) 



where: 



(n' x n' y n' z \V\n x n y n^j 

= i {ny ~ n ' y) HI 0< {x)<t>n> y {v)<t>n' z (z)V(x, y, z)(f) nx (x)(j) ny (y)<p nz (z)dxdydz 



-oo +00 +00 

■(ny-n') 



jjj e-^y 2 +^h n ,(x).h n ,(y).h n ,(z) 



— OO —OO —OO 



x y z 

x V< ^Y ) W ^' hnx ^' hnv<KV ^' hn ^' dxdydz ^ 



with f3 x , fl y , /5 z defined in eq.(24) 

• Due to the parity of V , the integral (30) vanishes if one of the three following 
conditions is not fulfilled: 

(-1)™- = (-!)< 

\-\) n v = {-l) n 'y 
(-1)«« = {-!)< 

Therefore, it is not necessary to calculate them. 

Since the respective indices must have the same parity, the complex factor 
of (30) can be rewritten as: 

jin' y -n y ) _ {_l\(n y -n y )/2 

• Although the (n' x n' y n' z \V\ n x n y n^j elements are spin independent, they are 

stored actually as:(n x n' y n' z Ti' V \n x n y n z T) in the computer memory 

• The matrix elements of the Coulomb energy are calculated in the same way 
as those of V(~r) (putting e(jf(~r) instead V(~r')), with the same change of 
scale. 

• In the gaussian integration, we have to calculate the Hermite polynomials, 
the central mean potential V(r), the spin-orbit mean potential S(r), and 
the Coulomb potential cj) c (~r') only at the nodes. Therefore, it is more con- 
venient to store these quantities in specific arrays before any calculations 
(see the common/tabh/... declaration in the subroutine setsub). 
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5.3.2 Matrix elements of the spin-orbit energy operator V so 



Due to the presence of the derivative of SCr'), the direct calculation of the ma- 
trix elements of V so , i.e. (n' x n' y n' z T! V so \n x n y n z T) , is not convenient. These 
derivatives can be transferred on the basis functions by partial integration. 
Therefore, the derivatives of the basis functions can be expressed from the basis 
functions themselves by mean of the recursion relations. Finally, the matrix el- 
ements of V so can be obtained by suitable combinations of (n' x n' y n' z \ S\ n x n y n z ^j 
elements, i.e. 



rL x n! y ri! z Y! \ V so \ n x n y n z Il 

^ [2B Z (£' \a z \ £) + B + (£' \a_\ E) + 5_ (E' \a + \ E)] (31) 



where: 



a± = a x ± o y (32) 



B± = B x =f (33) 



^ = l ^±{-^ ny ( nz + l) _ l ; n ^ j^l fi.nj,^ + l) 



^6 



- \Jn y (n' z + 1) (n' x n' y n' z + 1 j^l n x rij, - 1, n z ) 
+ \Jn z (n' y + 1) + 1, n' z l^l n x n y n z - l) 

+ \/ n 'z( n v + 1) - 1 1^1^% + (34) 



y = \j-^[-\f^Mx + 1) k,njn' 8 - 1 15| ra x + l,n y n z ) 

+ \Jn z (n' x + 1) (< + 1, |S| n.^n, - l) 
- \jn x (n' z + 1) (n x n y n' z + 1 l^l n x - 1, n 2 ) 
+ \jn' x {n z + 1) - l,n y n' z \S\n x n y n z + l)] (35) 
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B z = J^h/niK+lj (< - |5 , |n x 7i w + l,n z ) 

- \/n x (n' y + 1) + 1< n x - 1, n y n 2 ) 

+ + 1) (n' x + l ,n' y n' z \S\ n x n y - 1, n x ) 

+ ^(n* + 1) ( n X - 1, < 1^1 n x + 1, n v n x) (36) 

where S'(~r') is given by eq.(7) 

• The changes of sign in B y are involved by the phase factor of the basis 
functions. 

• The computations of the matrix elements of S(~r) are carried out like those 
of VCr)) in (30). 

5. 3. 3 Matrix elements of the kinetic energy operator T 

At last, the matrix elements of the T operator can be calculated in straight- 
forward way: 



(n' x n' y n' z Y,'\T\n x n y n z Y) = -cW^^x 

[f^ Z S nx n'Jn y n' y Sn z n' z (2n z + 1) 

- fku z 5 nxn i x 8 nyn > y 8 nzn > z+2 \Jn' z {n z + 1) 

- fkV z Sn x n'Jn y n' y S nz n' z -2 \J 'n z (n' z + 1) 

+ cyclic permutations] (37) 



6 Symmetry properties of the nuclear surface 

The three surfaces can be written as: 

iP' z^ 

7r = ^ + ^ + c^" 1 = (38) 

This implies the following properties: 

n(x, V, z) = ir(-x, y, z) = n(x, -y, z) = n(x, y, -z) (39) 
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Thus, for the two mean fields (i.e. central and spin-orbit fields), and for the 
Coulomb potential, we obtain: 

V(-x,-y,-z)=V(x,y,z) (40) 



S(-x,-y,-z)=S(x,y,z) (41) 



$ c (-x,-y,-z) =$ c (x,y,z) (42) 

6. 1 Parity 

Because of the relations (40), (41), and (42) the parity is a good quantum 
number, and the initial matrix decays into two sub-matrices according to the 
number 

p a = = ± i (43) 

Obviously, if (n x + n y + n z ) is even or odd the parity is respectively positive 
or negative. 

6.2 Signature 

Furthermore, the Kramers degeneracy is expressed here, by the fact that the 
eigenvalues are doubly degenerated relatively to the signature quantum num- 
ber q K , which is defined by: 

q K = (-l) n * +n yj: = ±1/2 (44) 

Consequently, the secular matrix splits into two sub-matrices, and only one 
must be considered. The two matrices contain the same set of eigenvalues, but 
the eigenf unctions are time-reversed each other. 

6.3 Consequences of these symmetries 

The computer code carries out calculations only for one kind of particles. 
Therefore, in order to take into account both neutrons and protons, the code 
must be run twice. 
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Since the Hamiltonian connects only states with the same parity, the computer 
code is built in such way that it separates the two types of parity p a — ±1 and 
performs the calculations separately for them. Consequently, the representa- 
tive matrix of the Hamiltonian splits into two blocks with a definite parity for 
each block. The diagonalization is then carried out in each block. 
Furthermore, the Kramers degeneracy involves the same eigenvalues for states 
which are time-reversed each other. For each block of a definite parity, the 
eigenenergies can be separated into two sets defined by the signature qx = 
±1/2. The code will make calculations only for q K = +1/2. The second block 
q K = —1/2 will be implicit, and will contain same energies but with time- 
reversed eigenfunctions. These eigenvectors may be obtained by application of 
the time reversal operator, i.e. by the operator T = —ia y Ko , where a y is a 
Pauli matrix and K is the operator of complex conjugation. 
Thus, one obtains 8 blocks, of which 4 are actually calculated (i.e. here the 
four first). 



1) 


M 


[Pa 


= +1] 


[Qk 


= +1/2] 


2) 


M 


[Pa 


= "1] 


[Qk 


= +1/2] 


3) 


[p] 


[Pa 


= +1] 


[Qk 


= +1/2] 


4) 


[p] 


[Pa 


= "1] 


[Qk 


= +1/2] 


5) 


[n] 


[Pa 


= +1] 


[Qk 


= -1/2] 


6) 


M 


[Pa 


= "1] 


[Qk 


= -1/2] 


7) 


[p] 


[Pa 


= +1] 


[Qk 


= -1/2] 


8) 


[p] 


[Pa 


= "1] 


[Qk 


= -1/2] 



So, it is important to point out that, the number of the basis states practically 
taken into account by the code is the half of the actual number. 



7 Numerical choices and prescriptions 

7.1 The quadratures 

The matrix elements of V( r ), e$ C7 ( r ) and V so ( r) are calculated with the 
Gauss-Hermite method, with 30 x 30 x 30 of mesh points. The Coulomb po- 
tential $ c ( ~r) is also evaluated numerically by the Gauss- Legendre method, 
but with 48 x 48 of mesh points . 

These choices seem to be sufficient relatively to the size of the basis ( iV max ^ 
26 ), and the interval of deformation ( < (3 < 0.6). A direct checking has 
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been done by increasing the number of quadrature points and by comparing 
the stability of the results (even with 20 points the results remain very correct). 



7.2 Prescription of the basis truncation 

In practice, the Hamiltonian matrix is finite. Therefore, for reasons of accuracy, 
we have to select a sufficient number of the basis states. Generally, we adopt 
one of the two following criteria: 

The first (spherical criterion) consists in choosing all basis states which satisfy 
the following inequality. 

n x + n y + n z < iV max (45) 

With this criterion the total number of the basis states is given by (Nmax+1) 
x (Nmax+2) x (Nmax+3) /6. 

In the second criterion (deformed criterion), one selects the states according 
to the deformation of the basis, i.e. according to the three frequencies of basis. 

Ill 3 
(n x + -)fkj x + (n y + -)fku y + (n z + -)ftw z < E cut = (N max + -)fku (46) 

(In fact, these three frequencies are already connected by the condition (27)). 
Thus, the choice of iV max determines the size of the basis. The files "con- 
vert. res" and "converl3.res" give some details about this. 



7. 3 Optimization of the basis frequencies 

Since the Hamiltonian operator does not depend on the oscillator frequencies, 
its eigenfunctions, and its eigenenergies, must not depend on these parameters. 
In practice, the representative matrix of the Hamiltonian is built by means of a 
finite number of oscillator eigenfunctions. This implies a spurious dependence 
according to these parameters. 

In another point of view, we might consider this method as a variational 
method in which the variational parameters are the frequencies of the basis. 
Thus, the best set (in terms of energy) for these frequencies should be precisely 
the one, which minimizes the eigenenergies, or simply their sum. 
For practical reasons, this method is not easy, since the variation is three- 
dimensional. However, it can be often more efficient to use some prescriptions 
in order to find (in an economical way) suitable values for these parameters. 
In the present work, we have adopted the approach of the references [1] and 
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[9]. In that method, we define first, the quantities p and q by: 

2 = = JdT.p(r)z 2 2 = (z 2 ) = JdT.p(r)z 2 

{x 2 ) JdT.p(r)x 2 (y 2 ) J dr.p(f)y 2 

where p(r) is the nuclear density. Note that the present definition of p differs 
from that of the ref.[l]. 

For a harmonic oscillator, the equations (47) are reduced to very simple 
relations. 

quo = — Pho = — (48) 

Next, we have to add, to these two formulas, the relation (27). Now, it is possi- 
ble to replace the parameter set ( u x , u y , u z ) by the equivalent (quo, Pho, ^o)- 
In the same way, for the potential of Woods-Saxon, the nuclear density can 
be approximated by the one of the liquid drop (i.e. a constant density). We 
obtain thus: 

c c 
qws = - Pws = 7 (49) 

a b 

At last, we "adapt" the oscillator basis to the nuclear shape by requiring: 

Qho = Qws Pho = Pws (50) 

For the u value, we can adopt simply the one of the Nilsson model. 

/lo^o ~ 41.A~^ (51) 

Many tests have shown that relations (50) and (51) give automatically very 
close values to those that produce the "true" minimization. Furthermore, a 
general rule is that a large basis size involves always a weak dependence of 
the eigenvalues according to these parameters. Going to the limit, we can say 
that if the basis was infinite, the results would be independent to the basis 
parameters. Conversely, for a too small basis, the dependance is strong, and 
the results become too inaccurate. 

For a square well, or (approximately) a Woods-Saxon potential, simple ana- 
lytical considerations lead to a more "refined" value for the parameter Jtlu : 



^o-t; ( — ) | Vb| -A' 1 / 3 w 0.979 \V \ .A~ 1/3 (52) 



5 (2\ 1 ' 3 
3 U 2 

where Vq is the depth of the potential. The equation (52) is obtained by 
requiring the condition 

(r 2 ) = (r 2 ) 

\ I harm. Oscillator \ / square well 

The averages are made with the semi-classical Thomas-Fermi density: 

"- w = w (A - l/(r))3/2 
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The Fermi level A is determined by the condition of conservation of the particle 
number A. 

The relation (52) could explain the empirical scale factor used sometimes 
[1,2,5] in the " standard equation" (51). 



7.4 The numerical values of the (3 parameters of the basis 



The quantities (3 = ^ , (3 X = ^ , (3 y = ^ , and (3 Z = ^p. are 
numerically calculated like 



muj / mcr 



The values are: 
m p c 2 = 938.2592 MeV 
m n c 2 = 939.553 MeV 
he = 197.32879 MeV fm 
This involves: 



m n c 



2 



" -t p = 0.0240958315 MeV^fm" 2 (53) 



(hey 



9 

777 C 

= t n = 0.0241290571 MeV' x fm~ 2 (54) 

(he) 2 



so that : 



P = y/t.Hvo, /3 X = y/t.hvx, f3 y = yjt.tkjy, (3 Z = \Jt.ftw z (55) 
with t — t p or t — t n 

The numbers t p and t n appear in the subroutine "Basisparam" . 



8 Diagonalization 



The diagonalization of the representative matrix of the Hamiltonian is carried 
out by a set of subroutines extracted from the EISPACK library of FORTRAN 
programs (http//www. netlib.org/eispack/). Thus, four subroutines of this li- 
brary were gathered: 

The subroutine tredl transforms any full symmetrical matrix into a tridiago- 
nal symmetrical matrix by using the Givens-Householder's method. 
For a tridiagonal symmetrical matrix the subroutine tqll uses the ql method 
to calculate only the eigenvalues of a tridiagonal matrix. 
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For a tridiagonal symmetrical matrix, the subroutine tsturm calculates the 
eigenvalues contained in a given interval. This subroutine calculates also the 
eigenvectors associated to the found eigenvalues. The adopted method is that 
of the bisection and the inverse iteration. 

Lastly, the subroutine trbakl recalculates the eigenvectors found by tsturm 
relatively to the initial basis (that of tredl). The sought eigenvectors are thus 
obtained. 

These subroutines are called by the subroutines diagoplus (for the positive 
parity), and diagominus (for the negative parity) in which the options of the 
diagonalization are specified. These options are indicated in the comments of 
the program, and below, in the subsection 10.1.1. 



9 The subroutines and the functions of the Program. 

The program is composed by a main program, 29 subroutines and 6 functions. 
The role reserved for each program is briefly described in the paragraph below 
(and described again in details in the comments of the program). 
In fact, all calculations are governed by the subroutine setsub which is in some 
sense a super subroutine. 

9.1 The set of subroutines (in the order of the calls) 

(1) The subroutine readl: reads the basic input parameters in the file in- 
put, dat 

(2) The subroutine writel : performs some tests and writes on a files eigvals.res 
and conver.res 

(3) The subroutine setsub: drives the successive calculations 

(4) The subroutine write2: writes on the file eigvals.res 

(5) The subroutine writeS: writes on the file conver.res 

(6) The subroutine woodsparam: calculates the Myers parameters 

(7) The subroutine surfparam: calculates the surfaces parameters 

(8) The subroutine basisparam: calculates the oscillator basis parameters 

(9) The subroutine pottablo : stores the potential at the nodes of quadrature 

(10) The subroutine coefftablo: stores the products of the coefficient of quadra- 
ture 

(11) The subroutine hermitablo: stores the Hermite polynomials at the nodes 

(12) The subroutine coultablo: stores the coulomb potential at the nodes 

(13) The subroutine statesplus: selects the numbers and the oscillator basis 
states corresponding to the positive parity 

(14) The subroutine statesminus: selects the numbers and the oscillator basis 
states corresponding to the negative parity 

(15) The subroutine idm: calculates the total numbers of the used basis states 
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(16) The subroutine matpotplus: calculates the representative matrix of the 
central mean potential for the positive parity 

(17) The subroutine matpotminus: calculates the representative matrix of the 
central mean potential for the negative parity 

(18) The subroutine matcinplus: calculates the representative matrix of the 
kinetic energy for the positive parity 

(19) The subroutine matcinminus: calculates the representative matrix of the 
kinetic energy for the negative parity 

(20) The subroutine matpotsoplus: calculates the representative matrix of the 
mean spin-orbit energy for the positive parity 

(21) The subroutine matpotsominus: calculates the representative matrix of 
the mean spin-orbit energy for the negative parity 

(22) The subroutine diagoplus: diagonalizes the representative matrix of the 
hamiltonian for the positive parity 

(23) The subroutine diagominus: diagonalizes the representative matrix of the 
hamiltonian for the negative parity 

(24) The subroutine tredl : Eispack subroutine (see section 8) 

(25) The subroutine tqll : Eispack subroutine (see section 8) 

(26) The subroutine tsturm: Eispack subroutine (see section 8) 

(27) The subrouine trbakl : Eispack subroutine (see section 8) 

(28) The subroutine eigenvalues: gathers the eigenvalues for both parities 

(29) the subroutine vektors: writes the eigenf unctions in a file. 



Fore several subroutines, the names ending in "plus" or in "minus" means 
that the subroutine performs calculations specifically for a defined parity. The 
term "plus" is employed for the positive parity, and the term "minus" for the 
negative parity. 



9.2 The set of functions 

(1) The function Hermite: calculates the Hermite polynomials 

(2) The function delta: delta symbol of Kroneker 

(3) The function potenv: calculates the central mean potential value at any 
point. 

(4) The function potenso: calculates the spin-orbit mean potential value at 
any point. 

(5) The function ephi: calculates the Coulomb energy of the proton at any 
point. 

(6) The function epslon: estimates the round-off error for the Eispack sub- 
routines 
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10 Input-output data of the FORTRAN program 

If no modifications are made the use of the program as presented in long 
theoretical description is very simple. 

10.1 The input data 

All input data are read from two files in a namelist type declarations. The 
second file is needed only if one does use a personal parameters for the po- 
tential, instead those of Myers. In this latter case, one has to precise its own 
parameters, in a second separate file. 

10.1.1 The first input data file: input.dat 

The file input.dat gathers all basic input data. Their significance is given 
below. 

• nmaxl and nmax2 are the bounds of the loop for iV max (eq(45) or (46)). This 
latter is the number of the major shells used in the calculations. If nmaxl 
=nmax2 (=nmax) the calculations are performed once. The variation of 
nmax is envisaged only if one desires to study the convergence of the results 
as a function of the number of the basis states. 

• pi is the pi number (3.1415927410125d.0) 

• If kkind—1, calculations are made for the neutrons case. 
If fcA;me?=i?,calculations are made for the protons case. 

Any other value of the kkind parameter involves an error declaration of the 
program. 

• Iz = number of protons. 

• In = number of neutrons. 

• Betta, and gama are the usual deformation parameters of Bohr (eq.(18a)- 
(18c)). 

• If ibase=0 the states of the basis are selected according to the spherical 
criterion (45). 

If ibase—1. The states of the base are selected according to the deformed 
criterion (46) There is not other value for this parameter. 

• If ili2=l, the program gives all eigenvalues, without eigenvectors. 

If ili2=2,the program gives the eigenvalues included in a given interval 
[elow, ehigh] with the corresponding (orthonormalized) eigenvectors. Any 
other value of this parameter involves an error declaration of the program. 

• Elow = lower bound of the selected interval. 

• Ehigh= higher bound of the selected interval. 

( Naturally, if this interval is sufficiently large it will contain all eigenvalues. 
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Consequently all eigenvectors will be also given.) 

• If icalc=0, the parameters of the Woods-Saxon potential are read from the 
namelist of the second input file parameters dat. 

if icalc—1, the Myers parameters are calculated by the subroutine woodsparam. 

• If iscal—1 the basis parameter Huq is computed from (52) i.e. from Huo = 
0.979 \V \. A' 1 / 3 

If iscal=2 the basis parameter Huo is computed from the relation Hujq = 
faktor.A^ 1 / 3 .(see eq.(51). 

• faktor= input parameter of the previous relation 

10.1.2 The second input data file: parameters.dat 

There is an option ( governed by the keyword icalc ) in the first input file 
which permits to the user to employ its own parameters instead of those of 
Myers. 

The data of the file parameters.dat are : 

• v0neut= deep of the central part of the potential for the neutrons 

• avneut= diffuseness of the central part of the potential for the neutrons 

• rvneut= radius of the central part of the potential for the neutrons 

• capasoneu= spin-orbit coupling strengh for the neutrons 

• assoneu= diffuseness of the spin-orbit part of the potential for the neutrons 

• rssoneu= radius of the spin-orbit part of the potential for the neutrons 

• v0pro= deep of the central part of the potential for the protons 

• avpro= diffuseness of the central part of the potential for the protons 

• rvpro= radius of the central part of the potential for the protons 

• capasopro= spin-orbit coupling strengh for the protons 

• assopro= diffuseness of the spin-orbit part of the potential for the protons 

• rssopro= radius of the spin-orbit part of the potential for the protons 

• rchpro= radius of the coulomb potential 

10.2 The output data 

The global results can be extracted from the five arrays evalplus, evalminus, 
evecplus, evecminus, and energies, in the main program. 
The arrays evalplus and evalminus contain respectively, the eigenvalues for 
the positive parity and the negative parity . The eigenvalues are classified in 
an increasing order. 

In the same way, the arrays evecplus and evecminus contain the components 
of the eigenvectors, in columns, in the same order as that of the eigenvalues. 
For the positive parity ( respectively the negative parity) , the parameter 
nevalplus in the subroutine diagoplus (respectively nevalminus in the subrou- 



21 



tine diagominus) gives the number of eigenvalues. 

Sometimes, it is more convenient to gather all eigenvalues in a common array 
(but the eigenvectors remain in their respective blocks). This is carried out in 
a common array named energies. In this array, the eigenvalues are classified 
in an increasing order. 

In order to find the corresponding eigenvector to a given eigenvalue, a vector 
containing a supplemental information was created and named num(k). The 
sign and the absolute value of num(k) indicate respectively the block (i.e. 
evecplus or evecminus) and the place of the column in this block. 

Furthermore, the output data can be consulted in a straightforward way, in 
three files: 

a) The eigenvalues are written in the file " eigvals.res ". In this file, it is 
indicated in particular, if the eigenvalues belong to the set corresponding to 
the positive parity or those corresponding of the negative parity. 

b) The eigenfunctions are recorded in the file " vekt.res ". For every eigenvalue 
there is a set of components relative to the different states (nx, ny, nz, sigma) 
of the basis. 

c) A brief study on the convergence is made in the file " conver.res " . 



11 Checking of the computer code and comments on the test run 

In order to check the code, one has proceeded to three types of tests. In the 
first, we use well-known analytical results. In the second, we compare our 
calculations with those using the same model. At last, in the third, we use 
some well-known properties of symmetry. 

11.1 Analytical tests 

In fact, the method of resolution of the Schrodinger equation proposed here is 
a purely numerical method. Consequently, one can use, not only the Woods- 
Saxon potential, but also any other type of potential. It is then possible to 
replace the Woods-Saxon potential by that of the harmonic oscillator in order 
to test the code by well-known analytical results. 

11.1.1 The deformed case without spin-orbit term 

Indeed, for a pure deformed harmonic oscillator, (without spin-orbit interac- 
tion), in cartesian coordinates, the theoretical expression of the energy is given 
simply by: 
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E(n x , n y , n z ) = (n x + \)fw x + (n y + \)tw y + (n z + \)tiw z 
n x ,n y ,n z = 0, 1,2, oo 

For reasons of simplicity, we have chosen the same frequencies as those of the 

basis. The numerical values are extracted from the file "eigvalsl.res" . 

In the computer program, one must replace the Woods-Saxon potential by 

that of the harmonic oscillator, i.e. by: 

V(~r) = \m(uj x x 2 + ujyi) 2 + oj 2 z z 2 ) 

Then, one has to cancel the spin-orbit interaction (by making the function 
potenso = or by cancelling the spin-orbit coupling constant) in the code. 
Calculating some levels analytically, and comparing them with those of the 
code, one can note an excellent agreement (to seven significant digits) for the 
deformed case (see table 1) 

It is important to note that the matrix elements are integrated numerically 
in the computer code, therefore, from this test, we can conclude that the pro- 
gram performs this task correctly. Because this term is diagonal ( in fact the 
code " does not know this" but after calculations, it finds that the nondiagonal 
elements are equal to zero ) in our basis, this test does not permit us to verify 
the diagonalization. These latter part of the program will be verified in the 
subsections below. 

Now, if we add the spin-orbit interaction, we could test the program entirely. 
Unfortunately, in the deformed case, there is no theoretical expression for that. 



11.1.2 The spherical case without, and, with spin-orbit term 



Of course, for the spherical case, it is possible to make fwo x = Hu> y = Two z = Huj 
in the previous theoretical expression. Nevertheless, the spherical coordinates 
are more convenient because as we shall see, the spin-orbit term has to be 
"treated" in that system. In this latter, the theoretical expression of the energy 
of a pure oscillator is well-known: 



2(n -!)+£ + 



D + 



fko = N+% hujQ 



This number specifies a major shell 



E(n, I, m) 
N = 2(n- 

n = 1,2, oo 

£ = 0,1,2, oo 

m = -1,-1 + 1,..., I, (for a fixed £) 

Due to the fact that the spherical symmetry is a particular case of the deformed 
case, it is obvious, that the results of the code (eivals2.res ) should be in 
complete agreement ( with the awaited degeneracy) with the analytical results. 
We can see in. table 2a,. that it is really the case. 

Furthermore, there, contrary to the deformed case, it is possible to obtain the 
analytical expression for the spin-orbit term. 

Indeed, one can use the relation (22) in a suitable way in order to obtain a 

simple theoretical expression for the spin-orbit term. 

Taking S(~r) = cr 2 /2, where c is a positive constant, one gets then: 



r or 



-c £ .a = -2c £ .~s = -2cUj 



2 ) 
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The most important point is that, in this way the splitting of the major shells 
does not depend on r, and, is then rigorously given by: 



AE{£-\) = c.{£+\) 
AE{£ + \) = -cl 



if 



if 

3 



3 



+ 



_ i 

2 

and 



and £ ^ 
£^ 



Thus, the new energies can be written as: 

E(n, I, j = £ T \) = [2(n - 1) + 1 + §] hu + AE(£ t \) 

Therefore, the code can be verified in its integrality. 

In order to simplify the numerical values of the splitting, we take c = lMeV. 
Thus, except for the value £ = 0, we can see that the levels are simply shifted 
by integer values according to the value of £. In order to illustrate that, we 
will give two examples: 



• Example 1 : 

if £ = 1, the p shell with a energy noted E(p) splits into two subshells 

according to the two values of j: 

for j = i + \ = 1 + \ = 3/2, 

E(p3/2) = E(p) -£ = E(p) - lMev 

forj=£-i = l-i = l/2, 

E(pl/2) = E(p) + {£ + 1) = E(p) + 2Mev 

• Example 2: 

similarly, if £ — 3 (/ shell ) one obtains 

for j = £ + \ = 3 + | = 7/2, 

E(f 7/2) = E(f) - £ = E(f) - 3Mev, 

for j = I - \ = 3 - \ = 5/2, 

E(f5/2) = E(f) + (£ + !) = E(f) + AMev 



In the table 2b, we compare all results of the code (eigvals3.res) with those 
of the analytical expression. Practically, the code ( which works in double 
precision) gives the exact values to six or seven significant digits for all levels 
of the spectrum. 

This high accuracy is due to the fact that the oscillator potential is a polynome 
of order two, therefore the Gauss method gives in this case the exact values 
for all matrix elements. 

Of course, these tests are not realistic cases, but they prove that the code runs 
properly with a high degree of precision. 

The case with spin-orbit term is very important because it involves the integral 
analytical checking of the code. Due to the fact that this operator is not 
diagonal in the oscillator basis, it proves not only that the code performs 
correctly all calculations of the matrix elements, but also proves that the step 
of the diagonalization is done properly. 

One also made some additional easy checks (not shown here). For example, by 
taking a constant potential in the spin-orbit term , one cancels the spin-orbit 
potential . That was well verified by the code, etc... 
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11.2 Comparisons with similar works 



For the deformed Woods-Saxon potential , it seemed to us more convenient 
to compare our code with those of the reference [2]. The reasons are the fol- 
lowing: 

a) We use exactly the same model as this reference. 

b) All potential parameters of the calculations are precised in that reference, 
and we need to use the same. 

c) Not only a part, but the entire spectrum of eigenvalues is given (as a func- 
tion of the deformation). 

The only disadvantage is that the results are displayed under a graphical form. 
However, in extracting the numerical values, we have tried to minimize the 
errors by using a graphical software. 

The eigenvalues are read with the own scale of the software. Then, a suitable 
linear transformation returns these values in MeV. Nevertheless, in order to 
find the "best values", this transformation has been carried out by the least- 
squares' method of the software. 

It turns out that it is possible to obtain values with an error about ±0.03MeV. 
We have thus considered the deformation ((3 = 0.3, 7 = 0.0) for the lead Pb208. 
For the basis parameters (Nmax and ftivo), we tried to use in calculations, the 
same, in order to obtain, as much as possible, close results. For Nmax, the ref- 
erence [2] indicates that the matrices corresponding to the two parities have 
a dimension of about 160 states. Consequently the fixed value for Nmax was 
certainly Nmax=10. However, the value of Hoj really used by the code is not 
given. This reference indicates only that, for the spherical case, the theoret- 
ical relation ftwo = 55. A^ 1 ^ 3 is better than the standard theoretical relation 
hcuo = 41. A~ x l 3 . Nevertheless, the reference [9], claims that a practical value 
of the order of 45-48MeV (instead 55. MeV) gives a somewhat better results 
that these theoretical relations. Since the codes of these two references have 
been compared, it is probable that a common practical value was fixed. We 
endeavored " to guess " this value. After many tests, It turned out that the 
value 47 MeV gave a good agreement 

Our calculations were carried out successively with Nmax = 10 (as the cited 
reference), and Nmax = 26. Indeed, this latter value insures that the levels 
are calculated with about three or four significant digits near the fermi level, 
and obviously, all the more for lower levels (see the file "converl3.res"). They 
are thus practically independent of the choice of the basis parameters 
In the tables (3a-3b,4a-4b) which have been deduced from the files "eigvals4.res" , 
"eigvals5.res", "eigvals6.res" , "eigvals7.res", we show respectively all bound 
levels of the Pb208 for four cases: 

a) neutrons-prolate shape ( r y=0° ) f3 = 0.3), Nmax=10 

b) protons-prolate shape ( / ~f=0°,(3 = 0.3), Nmax=10 

c) neutrons-prolate shape (j—0°,/3 = 0.3)Nmax=26 

d) protons-prolate shape (^=0°,f3 = 0.3)Nmax=26 
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The levels were separated in two distinct blocks according to their parity. 
Of course, for a finite potential, the discrete positive energy levels do not rep- 
resent, a valid solution of the continuum (see ref.[10]), therefore, we shall drop 
them . 

In all cases, we can note that the energy levels are practically the same ones 
for the low part of the spectrum, but relative small differences appear in the 
upper part of the spectrum. 

These differences are more prounonced for Nmax = 26 that for Nmax = 10. 
The analyse of these results leads to the following conclusions: 

• The lowest levels of the spectrum converge systematically more quickly than 
the others. As one goes up in the spectrum the convergence is in general 
slowest, but there can be some rare exceptions. 

• A rapid convergence involves a weak dependence relatively to the basis 
parameters. The highest levels of the spectrum are thus more sensitive to 
the basis parameters. One can affirm that if the basis parameters of our 
code are close to those of the reference [2], they are not rigorously the same 
ones. 

• In fact, one noted that this remark is general. Indeed, a modification of any 
parameter (for example those of the potential) in the calculations produces 
a modification relatively more significant for the highest levels than for the 
lowest levels. For example if the radius of the mean potential (spherical case 
for the neutrons) varies from 7.36fm to 7.40fm (all other parameters being 
constant), the first level, and the Fermi level undergo variations of 0.05MeV, 
and 0.28MeV respectively. The "general rule" is thus that the lowest levels 
are most "stable". 

• Owing to the fact that we employ very similar parameters, our results with 
Nmax = 10 are "artificially" very close to those of the reference [2] (the 
mean deviations are about 0.05 MeV for all cases). Thus, our purpose which 
was to recover the same results is now reached. But, the word "artificially" 
means that for this small basis, both results are not enough accurate, al- 
though they are the same 

Indeed, it is clear that they will be actually less precise than those obtained 
with Nmax = 26. Significant differences appear in the top of the spectrum. 
In the file "converl3.res" , one can note that the Fermi level is stabilized to 
about 0.01 0.03 MeV only starting from Nmax — 15 ^ 16. Therefore, 
calculations with Nmax — 10 ^ 14 produce mediocre results. 

The rapid convergence of the lowest states is due mainly to the fact that the 
corresponding wave functions are very similar to those of the oscillator. This 
is not the case for the highest states where the wave functions are strongly 
oscillating, and where the edge effect of the potential is "felt". 
This can be easily noticed in the components of the eigenvectors, in the file 
"vektl4.res" . For example, concerning the first eigenvalue, only the compo- 



26 



nents corresponding to the lowest quantum numbers are important (see the 
components numbered 1, 2, 12, and, 59). 



11.3 Tests using some properties of the parametrisation 7) 

Two simple tests can be carried out to check the consistence of the program: 

In the first, one compares the spectra obtained with the deformations(/3, 7) 
and ((3, —7) This operation is in fact nothing other that a simple permutation 
of the axes 1 and 2 of the ellipsoid. Of course the two shapes are the same, 
consequently, the respective spectra must be identical. 

In the files "eigvals8.res" and "eigvals9.res" one can easily check that is re- 
ally the case with an astonishing precision. In particular, one can note in 
these files the permutation of values of the parameters t\w x (hbaromegx) , and 
huiy (hbaromegy) . 

In the second, one compares the spectra obtained with the deformations 
(f3, 7 = 60°) and (— /3, 7 = 0°). There also, this operation is simply a cyclic 
permutation of the three axes of the ellipsoid. Therefore, the spectrum must 
also remain unchanged. 

As for the previous case, this can be easily verified in the files "eigvalslO.res 
and eigvalsll.res". 



11.4 Tests of convergence 

In the files "converl2.res" and "converl3.res" , we have shown the convergence 
of the sum (of the single particle energy) of the first 126 neutrons levels of 
Pb208 for two deformations. The potential's parameters are those of the ref- 
erence [2]. 

This sum has converged to less than 1 Mev only starting from the values 
Nmax = 14 and Nmax = 16, respectively for the spherical and the deformed 
cases. 

This implies for the Fermi level, a convergence to 0.02 MEV and 0.01 MEV 
respectively for these two cases. However, for Nmax ^ 16, theses deviations 
depend still of the value Huo- Obviously, for higher bound states , the precision 
will be less. 

Everything depends on what one wants to make. So, for example , for the 
Strutinsky's shells corrrection the previous values seem to be sufficient. 
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Always concerning the Fermi's level (converl3.res), one notices in general that 
it increases in absolute value as Nmax increases, but sometimes, it happens 
that it decreases slightly (in absolute value). For example, in the spherical 
case, it passes from 8.510 to 8.502 when Nmax passes from 13 to 14. We 
can easily see that the dimensions of subspaces corresponding to the positive 
and negative parity do not vary simultaneously when Nmax varies by one 
unit. For example, when one passes from Nmax=13 to Nmax=14 only one 
subspace, namely the one with a positive parity, undergoes changes from 252 
basis states to 372 basis states. The other remains the same with 308 basis 
states. 

In our example, the Fermi level belongs to the subspace of negative parity, 
therefore, apparently, it should not have to change. In fact, the formulae of 
the spin-orbit interaction connects the matrix elements of the two subspaces 
(see eq.34-36). This implies always a slight modifications in the subspace which 
has not varied, and this must not be assimilated to a noise. 

In fact, in this method, the "true noise" has two main sources : 

a) under-estimations of the number of points in the numerical integrations of 
the matrix elements. 

b) a too small basis or really inadequate values of the basis parameters. 

With 30 points of quadrature, a double precision, and a large basis (Nmax 
up to 26) these two problems are here minimized. 



12 Conclusion 



We have elaborated and checked a calculation program solving the equation 
of Schrodinger for a deformed potential of Woods-Saxon type. 
The program appears very rapid, and consequently, it becomes possible to use 
significant basis sizes. 

Calculations with small bases, like those which were carried out in the past 
with Nmax = 10 ^ 12 lead to a very poor precision. Our conclusions are 
corroborated by other works. For example, the ref. [5] has shown for Hartree- 
Fock calculations that the error in the energy of Pb208 is smaller than lMeV 
only for Nmax ^ 16. Other examples are given in the ref. [4] which confirm 
this fact. 

Similar codes [6,7] were made in the past, but with the assumption of axial 
symmetry. To our knowledge, triaxial Woods-Saxon calculations were never 
really undertaken with significant sizes of the oscillator basis. 
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A The Myers parameters 



The diffuseness parameters ay, a so , and the spin-orbit coupling k are the same 
as those of the Ref.[l]. 

a v = 0.66 fm (A.l) 
a so = 0.55 fm (A.2) 

k = 12 MeV fm 2 (A.3) 

The parameters of central potential, and of the spin-orbit potential, were ex- 
tracted from the Myers droplet model [3,8]. This theory uses Thomas Fermi's 
approximation to approach average properties of finite nuclei like the density 
radii, skin-thicknesses, in terms of neutron and proton numbers. 
In this model, two auxiliary quantities are first defined: 

3= ^ + aoii2^ (A4) 

1 + AW 

0.147 ^2 0.00248Z 2 



A l/S + - 3305 + ' A 4/3 ( A - 5 ) 

The physical significance of these quantities is explained in the Ref. [3] 
With help of these quantities, the depth of the mean potentials are written 
as: 

V (protons) = -52.5 - 48.75 (A.6) 

V (neutrons) = -52.5 + 48.75 (A.7) 

The radii of the central potentials (which are different for protons and neu- 
trons) are expressed by means of the nuclear density radii Ro(protons), or, 
Ro(neutrons), and the density diffuseness ay: 

R v (protons) = R (protons) { 1 - ^— ( - av ) \ (A.8) 

3 \R (protons) J ■ 

Rvineutrons) = R (neutrons) {l ( — — | )■ (A.9) 

' 3 \Ro(neutrons) J 1 



with 



56 ~ 

R (protons) = R + 0.82 - — - + 0.225 (A.10) 

Ro 

56 ^ 

Roineutrons) = R + 0.82 - — 0.225 (A.ll) 

Ro 
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R = r Q A l '\l-t) 



(A.12) 



r = 1.16 fm 

The radius of the spin-orbits mean field is given in the same way: 



(A.13) 



(A.14) 



At last, the radius of the charge density is given by: 




1/3 



( 



N- Z 



(A.15) 



N + Z 
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